MATLAB求解矩阵特征值的六种方法

2024-06-25 14:05| 来源: 网络整理| 查看: 265

MATLAB求解矩阵特征值的六种方法

关于这个特征值的求解一共六种方法幂法反幂法 QR方法对称QR方法 jacobi方法二分法

接下来就着重讲解这些算法的是如何使用的幂法算法如下，输入: 矩阵A、非零矢量x0、maxit(2000）、tol（1.0e-7) 输出：模的最大特征量a、模的最大特征量对应的特征向量x

function [a,x,n] = pmethod(A,x0,maxit,tol) if nargin == 3 tol = 1.0e-6; elseif nargin == 2 maxit = 1000; tol = 1.0e-6; end a0 = 0; y = A * x0; a = max(abs(y)); x = y / a; n = 1; while (abs(a - a0) > tol) && (n maxit) disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']); else disp(['幂法迭代',num2str(n),'步收敛']); end

反幂法应用幂法于矩阵的反求矩阵的特征值输入：矩阵A、非零矢量x0、maxit(2000）、tol（1.0e-7) 输出：模的最小特征量a、模的最小特征量对应的特征向量x

function [a,x,n] = vpmethod(A,x0,maxit,tol) if nargin == 3 tol = 1.0e-6; elseif nargin == 2 maxit = 1000; tol = 1.0e-6; end [L,U,P] = lu(A); a0 = 0; y = utri(U,ltri(L,P*x0)); a = max(abs(y)); x = y / a; n = 1; while (abs(a - a0) > tol) && (n = maxit) disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']); else disp(['反幂法迭代',num2str(n),'步收敛']); end

这里还用到上次的ltri和utri函数，使用时把他们当作函数来使用就可以

QR方法利用正交相似变换输入: 矩阵A、maxit(2000）、tol（1.0e-7) 输出：所有特征值、所有特征值所对应的特征向量

function [a,x] = qrmd(A,maxit,tol) %a(i,1)为第i个特征值，x(:,i)为第i个特征值对应的特征向量 if nargin == 2 tol = 1.0e-6; elseif nargin == 1 maxit = 1000; tol = 1.0e-6; end A0 = A; a0 =diag(A); [Q,R] = qr(A); A = R*Q; a = diag(A); n = 1; while (max(abs(a-a0)) > tol) && (n maxit) disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']); else disp(['QR方法迭代',num2str(n),'步收敛']); n = size(a,1); x = zeros(n); for i = 1:n x0 = ones(n,1); [b,x0] = vpmethod(A0-a(i,1)*eye(n),x0,maxit,tol); x(:,i) = x0; a(i,:) = a(i,:) + b; end end

这里也是运用到了一些常用的计算函数，直接复制在上面就可以

%ltri function y = ltri(L,b) n=size(b,1); y=zeros(n,1); for j = 1:n-1 y(j) = b(j)/L(j,j); b(j+1:n) = b(j+1:n) - y(j) * L(j+1:n,j); end y(n) = b(n)/L(n,n); %utri function x = utri(U,y) n=size(y,1); x=zeros(n,1); for j = n:-1:2 x(j) = y(j)/U(j,j); y(1:j-1) = y(1:j-1) - x(j) * U(1:j-1,j); end x(1) = y(1)/U(1,1); %vpmethod function [a,x,n] = vpmethod(A,x0,maxit,tol) if nargin == 3 tol = 1.0e-6; elseif nargin == 2 maxit = 1000; tol = 1.0e-6; end [L,U,P] = lu(A); a0 = 0; y = utri(U,ltri(L,P*x0)); a = max(abs(y)); x = y / a; n = 1; while (abs(a - a0) > tol) && (n = maxit) disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']); else disp(['反幂法迭代',num2str(n),'步收敛']); end

对称QR方法输入：矩阵A、maxit(2000）、tol（1.0e-7) 输出：所有特征值、所有特征值所对应的特征向量

function [a,x] = wilkinsonqr(A,maxit,tol) %a(i,1)为第i个特征值，x(:,i)为第i个特征值对应的特征向量,A为实对称矩阵 if nargin == 2 tol = 1.0e-6; elseif nargin == 1 maxit = 1000; tol = 1.0e-6; end n = size(A,1); A0 = A; A = hess(A); a0 =diag(A); d = (A(n-1,n-1) - A(n,n)) / 2; u = A(n,n) + min(abs([d + sqrt(d^2 + A(n,n-1)^2),d - sqrt(d^2 + A(n,n-1)^2)])); [Q,R] = qr(A - u * eye(n)); A = R*Q + u * eye(n); a = diag(A); m = 1; while (max(abs(a-a0)) > tol) && (m = maxit) disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']); else disp(['隐式对称QR方法迭代',num2str(m),'步收敛']); n = size(a,1); x = zeros(n); for i = 1:n x0 = ones(n,1); [b,x0] = vpmethod(A0-a(i,1)*eye(n),x0,maxit,tol); x(:,i) = x0; end

这里是要用到已经学过的ltir、utri、 vpmethod函数，直接复制上面就好 jacobi方法最古老的方法输入：矩阵A、maxit(2000）、tol（1.0e-7) 输出：模的最大特征量a、模的最大特征量对应的特征向量x

function [a,V] = jacobi_eig(A,maxit,tol) %a的元素为A的特征值，V的第j列为对应于特征值a(j,1)的特征向量 if nargin == 2 tol = 1.0e-6; elseif nargin == 1 maxit = 1000; tol = 1.0e-6; end A0 = A; n = size(A,1); V = eye(n); b = A0 - diag(diag(A0)); w = sqrt(sum(sum(b .* b))); [r,c] = find(abs(b) > w); k = 1; while (w > tol) && (k c(i,1)) if (abs(A0(r(i,1),c(i,1)))

【本文地址】

公司简介

联系我们